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ABSTRACT 

We present two-dimensional line-of-sight stellar kinematics of the lens galaxy in the Einstein Cross, 
obtained with the GEMINI 8m telescope, using the GMOS integral-field spectrograph. The stellar 
kinematics extent to a radius of 4" (with 0.2" spaxels), covering about two-thirds of the effective (or 
half-light) radius Re — 6" of this early-type spiral galaxy at redshift z/ ~ 0.04, of which the bulge 
is lensing a background quasar at redshift Zg ~ 1.7. The velocity map shows regular rotation up to 
~ lOOkms^^ around the minor axis of the bulge, consistent with axisymmetry. The velocity dispersion 
map shows a weak gradient increasing towards a central value of ctq — 166 ± 2kms~^. We deproject 
the observed surface brightness from HST imaging to obtain a realistic luminosity density of the lens 
galaxy, which in turn is used to build axisymmetric dynamical models that fit the observed kinematic 
maps. We also construct a gravitational lens model that accurately fits the positions and relative 
fluxes of the four quasar images. We combine these independent constraints from stellar dynamics 
and gravitational lensing to study the total mass distribution in the inner parts of the lens galaxy. 

We find that the resulting luminous and total mass distribution are nearly identical around the 
Einstein radius Re = 0.89", with a slope that is close to isothermal, but which becomes shallower 
towards the center if indeed mass follows light. The dynamical model fits to the observed kinematic 
maps result in a total mass-to-light ratio Tdyn = 3.7 ± 0.5 Tqj (in the /-band). This is consistent 
with the Einstein mass Me = 1-54 x 10^" M© divided by the (projected) luminosity within Re, which 
yields a total mass-to-light ratio of T^ = 3.4 Tqj, with an error of at most a few per cent. We 
estimate from stellar populations model fits to colors of the lens galaxy a stellar mass-to-light ratio 
T^ from 2.8 to 4.1 Tqj. Although a constant dark matter fraction of 20 per cent is not excluded, 
dark matter may play no significant role in the bulge of this ~L* early-type spiral galaxy. 
Subject headings: gravitational lensing — stellar dynamics — galaxies: photometry — galaxies: 
kinematics and dynamics — galaxies: structure 
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1. INTRODUCTION 

In the cold dark m atter (CDM) paradigm for 
galax y formation (e.g. iKauffmann fc van den BoschI 
|2002|) . galaxies are embedded in extended dark mat- 
ter distributions wit h a specific and (nearly) univer- 
sal ra dial profile (e.g. iNavarro et al.lfl997l : iMerritt et al.l 
|2005() . Measurements of rotation curves from neutral 
hydrogen (H I) observations in the outer parts of late- 
type galaxies have provided evidence for the presence 
of dark matter in th ese systems already rn ore than 
two decades ago (e.g. Ivan Albada et al.l Il985f ). In the 
outer parts of early- type galaxies, how e ver, cold gas is 
scarc e (but see e.g. iFrairx et al.l Il994t iMorganti et al.l 
Il997t IWeiimans et al.l l2008f) . so evidence for a dark 
matter halo has to come from other tracers, such as 
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kinematics of sta rs, planetary n ebulae and globular 
clusters (e.g. iCarol lo et al.l 119951: iGerhard et al.l 120011: 
iRomanowskv et al.li 2003: C ote et al.ll2003[ ) or hot X-ray 
gas (e.g. iFukazawa et al.i i2006[ ). However, these trac- 
ers are not always (sufficiently) available, the observa- 
tions are often challenging, and the interpretation mod- 
eling dependent. Also in the inner parts of galaxies the 
amount, sh ape and profile of dark matter is still poorly 
known (e.g. lPrimackll2004 ). even though stellar kinemat- 
ics are readily available. 

A fundamental problem in using stellar kinematics (as 
well as other collisionless kinematic tracers) for this pur- 
pose is the mass-anisotropy degeneracy: a change in the 
measured line-of-sight velocity dispersion can be due to 
a change in total mass, but also due to a change in veloc- 
ity anisotropy. Both effects can be disentangled by mea- 
suring also the higher-order velocity moments (Deionghj 
fl987l : Ivan der Marel fc FranxlfT993t JGerhardl [l993l ^ . but 



only the inner parts of nearby galaxies are bright enough 
to obtain th e required high-quali t y kinematic measure- 
ment s (e.g. Iva n der Marel' 'l99ll IGerhard etlHI l2001l : 
Capp ellari et al . 2006). A unique alternative to break 
the mass-anisotropy degeneracy is provided by gravita- 
tional lensing. In case of strong lensing, the mass of a 
foreground galaxy bends the light of a distant bright ob- 
ject behind it, resulting in multiple images. From the 
separation and fluxes of the images the total mass dis- 
tribution of the lens galaxy can be inferred directly. The 
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velocity anisotropy then can be determined from the ob- 
served velocity dispersion, without the need for higher- 
o rder velocity moments. 

iTreu fc Koopmana ()2004l . and references therein) have 
applied this approa ch to several strong tensi ng systems, 
of which 0047-281 (jKoopmans fc Treull2003[ l is the best 
constrained case, with three radially separated velocity 
dispersion measurements extending to about the effec- 
tive radius of the lens galaxy. They measure the mass 
within the Einstein radius by fitting a singular isother- 
mal ellipsoid to the positions of the quasar images. This 
Einstein mass is used to set the amplitude of the to- 
tal (stellar and dark) matter distribution, which they 
assume to be spherical. The constant stellar mass-to- 
light ratio T^, determines the contribution of the stars 
with a fixed radial profile, with the remainder due to 
dark matter with a single power-law profile with slope 7. 
They then compare the dispersion profile predicted by 
the spherical Jeans equations, for an ad-hoc assumption 
of the velocity anisotropy /3, with the observed disper- 
sion measurements. Based on a reasonably constrained 
Ti, and an upper limit on 7, they conclude that a signifi- 
cant amount of dark matter is present in the inner parts 
of the lens galaxy, with a slope shallower than the nearly 
isothermal total mass distribution. An additional (exter- 
nal) constraint on T^, (or 7) is needed to go beyond this 
limit on the dark matter distribution. Even so, their re- 
sults are limited by too few kinematic constraints (which 
leaves the anisotropy degenerate), and by the use of a 
simple spherical dynamical model. 

Clearly, most lens galaxies are significantly flattened 
and so cannot be well-described by spherical models. 
Non-spherical models provide a more realistic descrip- 
tion of the lens galaxy, but the increase in freedom re- 
quires also (significantly) more spatially resolved kine- 
matic measurements to constrain them. Only a few of the 
known strong gravitational lens systems are close enough 
to obtain such kinematic data. Even then, one can make 
(ad-hoc) assumptions on the velocity distribution, such 
as velocity isotropy in the meridional plane of an axisym- 
metric model. The latter restriction to a distribution 
function of two (instead of three) integrals of motion al- 
lows for the recovery of the intrinsic shape and mass dis- 
tribut ion in the presence of two-d i mensional kinemati c 
data (jBarnabe fc Koopmansl l2007t |Czoskc et ar| l2008D . 
possibly even without the ad ditional c onstrain t prov ided 
by gravitational lensing (e.g. lCappellari et al.ll2006l ). 

In this paper, we relax the two-integral assumption 
and find that one can still derive independent determina- 
tions of the intrinsic mass distribution from both gravita- 
tional lensing and stellar dynamics with two-dimensional 
kinematic data. We have observed the gravitational 
lens system QSO 2237-1-0305, well-known as the Einstein 
Cross, with the integral-field spectrograph CMOS on the 
GEMINI-North Telescope. We fit the extracted steflar 
velocity and velocity dispersion maps of the inner parts 
of the lens galaxy at a redshift zi ~ 0.04 with axisym- 
metric dynamical models. The resulting intrinsic mass 
distribution agrees well with that inferred from the lens 
model that fits the four quasar image positions and rel- 
ative fluxes. Furthermore, we compare this total mass 
distribution with the stellar mass distribution, based on 
an independent estimate of the stellar mass-to-light ratio 
T^,, to place constraints on the dark matter distribution 



in the inner parts of the lens galaxy. 

In Section[2l we briefly describe the Einstein Cross, and 
we present the photometric and spectroscopic observa- 
tions. In Section [3l we extract the two-dimensional stel- 
lar kinematics, describe the dynamical modeling method, 
and construct lens models. From the latter we infer the 
total mass distribution, which we compare in Section [5] 
with the luminous mass distribution inferred from the ob- 
served surface brightness. We then build axisymmetric 
dynamical models, and compare the resulting dynami- 
cal mass-to-light ratio with that inferred from the lens 
models, and in turn with an estimate of the stellar mass- 
to-light ratio. In Section [5] we discuss our findings and 
we summarize our conclusions. 

Throughout we adopt the WMAP cosmological pa- 
rameters for the Hubble constant, the matter den- 
sity and the cosmological constant, of respectively 
Ho = 73kms-iMp c-\ 17m = 0.24 and ^l = 0.76 
(iSpers-elet al.l[2?)0l . 

2. OBSERVATIONS AND DATA REDUCTION 

2.1. The Einstein Cross 

The Einstein Cross is the well-known gravitational 
lens system QSO 2237-1-0305 (RA: 22h 40m 30.3s, Dec: 
-K03°21'31", J2000). In this system, a distant quasar 
at redshift Zg = 1.695 is lensed by the bulge of the early- 
type spiral PGC 069457 at zi = 0.0394 (angular diameter 
distance Di = 155 Mpc, 1" = 0.75 kpc), resulting in a 
cross of four bright images separated by about 1.8". 

The Einstein Cross has long been the closest strong 
gravitational lens system kn own, and has been v ery well 
studied since its discovery bv lHuchra et al.l (|1985f ). There 
is a wealth of grou nd- and space-base d imaging data at all 
wavelengths (e.g . iFalco e t al. 1996; Blanton et al.lll998t 
lAgol et al.ll2000l : lDai et''al.ll2003l) . The resulting precise 
measurements of the positions and relative fluxes of the 
quasar images can be used to construct a detailed lens 
model. 

In contrast, kinematic data of the lens galaxy is 
very scarce, with onl y one measured central stellar ve- 
locity dispersion (Foltz et al. 1992) and two H I ro- 
t at ion curve measurements in the very outer parts 
(jBarnes et al.llT999 ). There are several pr evious integral- 
field s tudies of the E instein Cross: TI G ER: iFitte fc AdamI 
(I1994D: INT EG RAL: iMediavilla et all ([19981); CIRPASS: 
iMetcalf et al.l ()2004f) . However, none of these studies 
were concerned with the stellar kinematics of the lens 
galaxy, but instead investigated the quasar spectra. 

2.2. Imaging 

We use two WFPC2 images retrieved from the HST- 
archive (F555W V-band image, 1600 seconds, PI: West- 
phal, and F814W /-band image, 120 seconds, PI: 
Kochanek; see left and middle panels of Fig. [T|) to 
determine the luminosity density of the lens galaxy. 
We correct the J-b and image for extinction following 
ISchlegel et al.l |1998") , and we conve rt to solar units using 
the WFPC2 calibration of lDolphinI HOOO), while assum- 
ing an absolute /-band magnitude for the Sun of 4.08 
mag (Table 2 of iBinnev fc Merrifieldl[l99l ). From a de 
Vaucouleurs R^^^ profile fit to the /-band photometry 
in the inner bulge-dominated region, we obtain an effec- 
tive radius Re — 6 ", which is c onsistent with previous 
measurements (e.g. lRacinill99lD . 
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Fig. 1. — Surface brightness as observed with HST (left and middle panel) and surface mass density of the overall best-fit lens model 
(right panel) of the lens galaxy in the Einstein Cross. Left panel: the contours of the WFPC2/F555W V-band image reveal clearly the 
bulge, spiral arms and bar embedded in the large-scale disk of this early-type spiral galaxy. The cUipticity measured from the MGE fit 
(contours) is used to estimate the inclination. Middle panel: the central 8"x8" of the WFPC2/F814W /-band image, of which the MGE 
fit is used to construct the (stellar) luminosity density model of the lens galaxy. We use the 7-band image instead of the longer exposed 
l^-band image as it tracers better the old stellar population and is less sensitive to extinction and reddening. The four quasar images are 
masked out during the MGE fit. In both WFPC2 images the contours are in steps of 0.5mag/arcsec^. The images are rotated such that 
North is up and East is to the left. Right panel: The scale-free lens model with slope a = 1.0 (dashed contours) fits the positions and 
relative fluxes of the quasar images, indicated by the filled circles. Superposed is the MGE fit (solid contours). 



For the construction of the lens model, we use the ac- 
curate positions of the quasar images from the website 
of the CASTLES survey^ based on Hubble Space Telescope 
(HST) imaging. Although also optical flux ratios are 
gi ven on this websit e, we use the radio fluxes provided 
bv lFalco et al.1 ()1996[ ). because they are in general (much) 
less affected by differential extinction or microlensing. 

2.3. Integral-field spectroscopy 

Observations of the Einstein Cross lens system were 
carried out usi ng the integral -field uni t of the GMOS-N orth 
spectrograph (jMurrav et al. 2003; H ook et al.1 I2004D on 
July 17*'' and August I''* 2005 as part of the program 
GN-2005A-DD-7. The data were obtained using the IFU 
two-slit mode that provides a field-of-view of 5"x7". An 
array of 1500 hexagonal lenslets, of which 500 are located 
1' away from the main field to be used for sky subtrac- 
tion, sets the 0'.'2 spatial sampling. Eight individual sci- 
ence exposures of 1895 seconds each were obtained during 
the two nights, resulting in a total on-source integration 
time of ~ 4 hours. An offset of O'.'S was introduced be- 
tween exposures to avoid bad CCD regions or lost fibers. 
The R400-G5305 grating in combination with the CaT- 
G0309 filter was used to cover a wavelength range of 
7800-9200 A with a FWHM spectral resolution of 2.8 A. 

For the data reduction we used an updated version of 
the officially distributed Gemini IRAF'^" package. We ap- 
plied bias subtraction, flat-fielding and cos mic ray rejec- 
tion (using the L.A. Cosmic algorithm bv Ivan DokkumI 
1200 It) to each science exposure, and wavelength calibra- 
tion after the extraction of the data. A careful flat- 
fielding procedure, crucial for proper removal of fringes in 
the spectral direction, was carried out using the Quartz 
halogen (QH) lamp. Wavelength calibration was per- 
formed with CuAr lamp exposures taken before each sci- 

® http://cfa-www.harvard.edu/castles/ 

^^ IRAF is distributed by the National Optical Astronomy Ob- 
servatories, which are operated by the Association of Universities 
for Research in Astronomy, Inc., under cooperative agreement with 
the National Science Foundation. 



ence exposure. At the observed wavelengths, a complex 
spectrum of H2O absorption features overlaps with the 
CaT lines we are interested in. We constructed a correc- 
tion spectrum from observations of the white dwarf star 
Wolfl346, taken with the same instrumental setup as our 
science frames. We checked the range of fiber-to-fiber 
variations of the spectral resolution of the instrument by 
measuring the width of the sky lines in each fiber, result- 
ing in the nominal FWHM value of 2.8 ± 0.2 A for most 
of the data cubes. In order to combine all the data, we 
homogenized the science frames by convolving all spec- 
tra to an instrumental FWHM resolution of 3.1 A (or 
fjinstr — 44.2 km s"'^) and resampling the spectra to the 
same range and sampling in wavelength. After interpo- 
lating each science frame to a common spatial grid, tak- 
ing into account the small spatial offsets applied during 
the observations, we sum the spectra sharing the same 
position in the sky to produce the final merged data cube. 
The significant contribution from the sky lines to the 
overall spectrum of the lens galaxy made the sky sub- 
traction the most challenging step in the data reduction 
process. Although the scatter in the instrumental res- 
olution was small within the individual science frames, 
the combined effect of other steps in the data reduction 
(e.g. fringing, telluric absorption, resampling of data in 
wavelength), introduced alterations in the shape of the 
sky lines (from the sky fibers) that were not equally re- 
produced in the science fibers. As a consequence any 
attempt to use the sky from the sky fibers resulted in 
serious systematic effects (i.e. P-Cygni like residuals) in 
the galaxy spectra underneath. In order to minimize 
these effects, we chose to extract the sky spectra from 
the outer most regions of our fully merged science data 
cube, where no galaxy light was appreciable. The sky 
was then subtracted optimally as described below. 

3. ANALYSIS 

3.1. Two-dimensional stellar kinematics 

To measure reliable stellar kinematics we first co-add 
spectra using the adaptive spatial two dimensional bin- 
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Fig. 2. — Spectrum from the centerer of the lens galaxy (black 
curve), showing the Ca II triplet region fitted by a composite of 
stellar population models (smooth gray curve). The peaks are the 
sky spectrum. 



ning scheme of ICappellari &: Copinl ()2003D to obtain in 
each resulting (Voronoi) bin a minimum signal-to-noise 
ratio (S/N). The resulting spectra are spectrally rebinned 
to steps of constant In A, equivalent to a sampling of 
24kms~^ per pixel. We adopt the single s tellar popu- 
lation (SSP) models of iVazdekis et al.l (|2003f ) as spectral 
templates, which are convolved to have the same effective 
instrumental broadening as the CMOS data. Next, using 
an updated version ( v4.2) of the penalized pixel- fitt ing 
(pPXF) algorithm of ICappellari fc EmsellemI ()2004D . a 
non-negative linear combination of these templates is 
convolved with a Gaussian line-of-sight velocity distri- 
bution (LOSVD) and fit to each observed spectrum, to 
derive the mean line-of-sight velocity V and velocity dis- 
persion a per bin on the plane of the sky. This version of 
the code allows for the inclusion of a sky spectrum into 
the set of stellar templates, and determines the optimal 
scaling of the sky to the overall spectrum. We used this 
procedure to clean the individual spectra in our merged 
data cube, before we derived the final stellar kinematics. 

As an example of how this routine performed, we show 
in Fig. [21 the spectrum from the center of the lens galaxy, 
with the best-fit SSP template overplotted. We also plot 
the removed sky spectrum on top, to show that the sky 
emission-lines can severely affect the Ca II triplet absorp- 
tion lines that we use to extract the stellar kinematics, 
especially at larger radii where the the relative contribu- 
tion from the sky is more significant. In particular, in 
our experiments, a poor sky subtraction caused signifi- 
cant systematic effects in the map of a at the receding 
side of the galaxy, where one of the absorption lines is 
redshifted onto (the residual of) a sky emission line. The 
map of y, however, was (nearly) unaffected, showing reg- 
ular rotation with a straight zero-velocity curve that co- 
incides with the minor axis of the (bulge) photometry. 
These tests already show that, at least in the inner parts 
of the lens galaxy, both the photometry and kinematics 
are consistent with axisymmetry. 

If we now assume this symmetry from the start and 
apply it during the extracting of the stellar kinematics, 
we can not only increase the S/N, but at the same time 
also suppress or even remove systematic effects. To do 
this, we make use of a built-in functionality in the pPXF 
algorithm that allows diametrically opposed spectra to 
be fitted simultaneously. The symmetry imposed by 
this technique can be either mirror-symmetric or point- 
symmetric, corresponding to the behavior of the LOSVD 



in axisymmetric or triaxial stellar systems, respectively. 
In practice, this i s don e following the method outlined 
bv iRix fc White! ()1992D . whereby the two spectra from 
opposed positions in the galaxy are laid end-to-end to 
create a single vector. The template library is also dou- 
bled in length, where the second half is convolved by the 
symmetric counterpart of the broadening function used 
for the first half (simply reflecting the broadening kernel 
around the systemic velocity). The two galaxy spectra 
are therefore fitted simultaneously, with the implied as- 
sumption of symmetry. 

Performing a symmetric extraction of the kinemat- 
ics presents a number of advantages, and has been 
used by numerous autho rs not just as a quality-check 
(e.g. iRix fc White! !l992f ). but also when the assump- 
tion of symmetry is a n integral part of the analysis (e.g. 
!Gebhardt et al.!r2003( ). The two main advantages for the 
case presented here are in reducing the impact of system- 
atic effects in the data (including sky subtraction) and 
increasing the effective spatial resolution. The latter is 
possible since under the assumption of axisymmetry, all 
moments of the LOSVD are symmetric about the galaxy 
major axis. We can therefore 'fold' the data cube about 
the major axis, combining spectra from the two hemi- 
spheres, which increases the effective S/N within a given 
spatial element, and therefore reduces the required spa- 
tial bin size. 

We interpolate the merged CMOS data cube onto a new 
spatial grid with x' and y' coordinates parallel to the 
major and minor photometric axes respectively. We then 
fold the data-cube about the major axis by summing to- 
gether spatial pixels ('spaxels') with equal (x, \y\) posi- 
tions. Before spatially binning the spectra to a minimum 
S/N, the folded data cube is folded again, into a single 
quadrant by combining spaxels with equal (|a;|,7;) posi- 
tions. The Voronoi bins are derived in this quadrant, and 
the resulting bin centers are mirrored around x' = 0. The 
spectra are combined according to these Voronoi bins to 
ensure symmetry around the minor axis. At the field 
perimeter, it is generally not the case that both sides of 
the galaxy are sampled at the exact same positions, due 
to the alignment of the CMOS field with respect to the 
galaxy's principle axes. This results in some bins where 
the S/N of the opposed spectra is far from equal, or in 
some cases where observations are only present on one 
side of the galaxy. The former is taken into account by 
consideration of the error spectrum of each bin. For the 
latter, the bin is fitted as a single spectrum, giving kine- 
matic values only at that position, without a mirrored 
counterpart. 

From Monte-Carlo simulations for data of this spectral 
range and resolution, a minimum S/N of 20 is deemed 
optimal. This results in an average error in both V and 
(J of about 9kms~^, while still preserving the spatial 
resolution of the data. However, towards the edge of the 
CMOS field, the errors increase to ~ ISkms""'^ in V and ^ 
20 km s~^ in a. The final symmetrically binned data cube 
consists of 4 X 71 = 284 bins, with the largest containing 
19 individual 0.2" spaxels. The resulting (original and 
symmetrized) maps of V and a are presented in Fig. [3l 
with superposed contours of the integrated flux from the 
data cube. The velocity map shows regular rotation with 
amplitude of T^ up to '^ lOOkms"^. The dispersion map 
shows somewhat higher values of a at the position of the 
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Fig. 3. — Mean line-of-sight velocity V and velocity dispersion cr maps of the lens galaxy in the Einstein Cross as measured from 
observations with the integral-field spectrograph GMDS on Gemini-North. At the right side of each map, the (linear) scale in kms^^ is 
indicated by the color bar, and the limits are given below. The third and fourth panel show the V and a maps after symmctrization 
assuming oblate axisymmetry (see ii l3.fl l. The overplotted (black) contours show the isophotes of the image reconstruction from collapsing 
the unsymmetrized data cube, highlighting the regions excluded in the symmetrized kinematics due to contamination from the quasar 
images. The V map shows clear and regular rotation around the (vertically aligned) short axis of the bulge. The a map shows a weak 
gradient decreasing towards the edge of the field. 



quasar images. This is caused by imperfect subtraction 
of the quasar continuum, which acts to strongly dilute 
the absorption lines and make them appear broader. 

In our stellar template fit to the observed spectra, we 
include additive polynomials to account for remaining 
systematic effects in the data, e.g. due to imperfect flat 
fielding. At the same time, we use these additive com- 
ponents to mimic the contaminating contribution of the 
quasar images, which impart a strong featureless con- 
tinuum on top of the galaxy's integrated stellar spec- 
trum. This method nicely removes nearly all contribu- 
tion from the quasar, except where the quasar images are 
the brightest, in which case a cannot be reliable mea- 
sured. We exclude these regions when we fit dynamical 
models. However, between the quasar images, the spec- 
tra are dominated by the stellar light of the galaxy, giving 
a reliable central dispersion (Tq — 166 ± 2kms^^. Like- 
wise, at the edges of the field, the quasar contamination 
is small, providing a reliable gradient in the velocity dis- 
persion. 

Comparing in Fig. [3] the symmetrized stellar kinemat- 
ics with what is obtained without imposing axisymme- 
try, there is a marked reduction of noise and systematic 
effects, especially in the velocity dispersion map. How- 
ever, the main features of regular rotation, decreasing 
a at larger radii, as well as corresponding velocity and 
dispersion amplitudes, are present in both approaches, 
albeit with larger uncertainty and systematics in the non- 
symmetrized case. Given the various advantages, we con- 
sider hereafter the symmetrically extracted kinematics. 

3.2. Dynamical models 

As mentioned above, both the photometry and (un- 
symmetrized) kinematics in the inner parts of the lens 
galaxy are consistent with axisymmetry. Hence, in con- 
structing dynamical models we consider an axisymmetric 
stellar system in which both the potential ^{R,z) and 
distribution function ( DF) ar e independent of azimuth 
(j) and time. By Jeans' (|1915f l theorem the DF only de- 
pends on the isolating integrals of motion: f{E,Lz,l3), 
with energy £■ = {v'j^ + v-j^ + vl)/2 + ^{R, z), angular mo- 
mentum Lz = Rv^ parallel to the symmetry z-axis, and 
a third integral I^, for which in general no explicit expres- 



sion is known. However, usually^^ I3 is invariant under 
the change {vii,Vz) —^ {—vr, —Vz). This implies that the 
mean velocity is in the azimuthal direction (Ur = TTJ = 0) 
and the velocity ellipsoid is aligned with the rotation di- 
r ection (vbVa = Vd >Vz = 0) . 

'Sch warzschildl ()1979f ) introduced a method that 
sidesteps our ignorance about the non-classical inte- 
grals of motion. It finds the set of weights of or- 
bits computed in an arbitrary gravitational potential 
that best reproduces all available photometric and kine- 
matic data at the same time. The method has proved 
to be powerful in building detailed spherica l and ax- 
isymnietric models of nearby galaxies fe.g. | Rix et all 
1997t Ivan der Marel et"al] 119981: iGebhardt et all 120031 



Wluri et al.ll2004l : ICappellari et al.ll2006l: iThomas et all 



2007) as well as globular clusters ( van de Ven et al.ll2006t 
van den Bosch et al. 2006), and since recently also triax- 
ial models ^van de Ven et al.1l2008at Ivan den Bosch etlHI 
.2008.) . However, in all cases the stellar systems are signif- 
icant closer than the lens galaxy in the Einstein Cross, al- 
lowing for (even) more and higher-order stellar kinematic 
measurements necessary to constrain the large freedom 
in this general modeling method. Instead we construct 
simpler, but still realistic dynamical models based on the 
solution of the axisymmetric Jeans equations. 

When we multiply the coUisionless Boltzmann equa- 
tion in cylindrical coordinates by respectively vr and Vz 
and integrate over all velocities, we obtain th e two Jeans 
equations (see also lBinnev fc Tremaindfr987l eq. 4-29) 



d{R>^vl: 



dR dz 

d{Ri^VRvl) 
dR 



R — ^ vvi 



5$ 



R 



dz 



Ru 



(9$ 

dz 



^0, 



(1) 

(2) 



where v{R^ z) is the intrinsic luminosity density. Due to 
the assumed axisymmetry, all terms in the third Jeans 
equation, that follows from multiplying by v^f,, vanish. 

We are thus left^ with four unknown second order ve- 
locity moments vj^, v^, v^ and vrv^ and only two equa- 
tions. This means we have to make assumptions about 

^^ If resonances are present, I3 may loose this symmetry. 
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the velocity anisotropy, or in other words the shape and 
ahgnment of the velocity ellipsoid. In case the velocity 
ellipsoid is aligned with the cyhndrical (i?, </), z) coordi- 
nate system vrVz = 0, so that we can readily solve equa- 
tion (m for u^ . If we next assume a constant flattening 
of the velocity ellipsoid in the meridional plane, we can 

write vfi = v'^/{l — f3) and solve equation H]) for vi. 
This assumption provides in general a good descriptio n 
for the kinematics of real disk galaxies (|Cappellarill2008f ) , 
such as the lens galaxy in the Einstein Cross, which is 
an early- type spiral galaxy. When /3 = 0, the velocity 
distribution is isotropic in the meridional plane, corre- 
sponding to t he well-known case of a two-in tegral DF 
fiE,Lz) (e.s. lLvnden-BeIlfl96l lHuntej[l97l . 

Knowing the intrinsic second-order velocity moments, 
the line-of-sight second-order velocity moment for a stel- 
lar system viewed at an inclination i > away from the 
z-axis follows as 



1 



^los- 



I{x',y')J- 



+ OC 



vl, SVC? I 



9 2 

vi, COS ( 



■ 2 • 

sm I 



' cos^ i 



vrVz sm ( 



)sin(2i) 



Az' 



(3) 



where I{x',y') is the (observed) surface brightness with 
the x'-axis along the projected major axis. For each po- 
sition (x', y') on the sky-plane, vf^^ yields a prediction of 
the combination V^ + g^ of the (observed) mean linc-of- 
sight velocity V and dispersion a. 

Under the above assumptions, besides the anisotropy 
parameter /3 (and possibly the inclination j), the only un- 
known quantity is the gravitational potential, which via 
Poisson's equation is related to the total mass density 
p{R,z). We may estimate the latter from the intrinsic 
luminosity density v{R, z), derived from deprojecting the 
observed surface brightness I{x',y'), once we know the 
total mass-to- light ratio Ttot- It is common in dynamical 
studies of the inner parts of galaxies to consider Ttot as 
an additional parameter and to assume its value to be 
const ant, i.e., mass follows light (e.g. CappcUar i et al.l 
|2006|) . Since Ttot may be larger than the stellar mass- 
to-light ratio T*, this still allows for possible dark matter 
contribution, but with constant fraction. As an alterna- 
tive to these 'standard' constant-Ttot models, we also 
construct dynamical models in which we use the strong 
gravitational lensing to constrain the gravitational po- 
tential. In S 13.31 below, we use the positions and relative 
fluxes of the quasar images to estimate the surface mass 
density T,{x',y'), which we deproject to obtain p(R,z), 
so that Ttot is not anymore a free parameter. For both 
type of dynamical models, the best-fits to the stellar kine- 
matic data are shown in Fig. [71 and discussed further in 
§|U below. 

We use the Multi-Gaussian Ex pansion method (MGE; 
iMonnet et al.lll99"3 lEmsellem et al. 1994) to parameter- 
ize both the observed surface brightness I{x', y') as well 
as the estimated surface mass density S(a;', y') as a set of 
two-dimensional Gaussians. Even though Gaussians do 
not form a complete set of functions, in general surface 
density distributions are accurately reproduced, includ- 
ing deviations from an elliptical distribution and elliptic- 
ity variations with radius. Representing also the point- 
spread function (PSF) by a sum of Gaussians, the con- 



volution with the PSF becomes straightforward. More- 
over, the MGE-parameterization has the advantage that 
the deprojection can be performed analytically once the 
viewing angle (s) are given. Also many intrinsic quanti- 
ties such as the potential can be calculated by means of 
simple one-dimensional integrals. Similarly, the calcula- 
tion of uj^^g in equation ([3]) reduces from the (numerical) 
evaluation of in general a triple integral to a straight- 
forward single integral (see CappcUari 2008, for further 
details and a comparison with more general dynamical 
models). The latter integral is given in Appendix \X\ to- 
gether with expressions for I'iR, z) and ^{R, z) for the 
MGE-parametrization of I{x',y') and E(x',y') as given 
in Table [2l and further discussed in § 14.11 below. 

3.3. Lens models 

Under the thin-lens approximation, the gravitational 
lensing properties of a galaxy are characterized by its 
potential projected along the line-of-sight, also known 
as the deflection potential (t){x\y'). In case of a MGE 
parametrization of the surface mass density of the galaxy, 
the projection along the line-of-sight of the correspond- 
ing gravitational potential, or, alternatively, solving the 
two-dimensional Poisson equation becomes straightfor- 
ward. The calculation of (j){x' ,y') and its derivatives re- 
duce to the (numerica l) ev aluation of a single integral, 
as shown in Appendix IA.2I for an axisymmetric system. 
This MGE approach thus provides a powerful technique 
to build general lens models. However, since in this case 
the quasar images only provide a few constraints within 
a limited radial range, we adopt instead a simpler and 
less general approach. 

We use the algorithm of lEvans fc WittI (|2003l ) to con- 
struct a lens model that accurately fits the (optical) po- 
sitions and relative (radio) fluxes of the four quasar im- 
ages in the Einstein Cross. The deflection potential is 
assumed to be a scale-free function (j){R, 9) ex R"F{9) of 
the polar coordinates R and 9 in the lens sky-plane, with 
< a < 2 for realistic models. The angular part F[6) is 
expanded as a Fourier series 



F{9) 



■Co 



E 

m—l 



c,n cos{m6) + Sm sin(m6')] . (4) 



The positions (xi,yi) = {Ri cos 9i,RiSm 9 i) of the images 
are related to the position (£, tj) of the source by the lens 
equation fe.g. [Schneider et al]|1992D 

^ = x,[l~ {aF - tan 9, F^) Kf" ^j , 

T^ = y, [I -{aF,- cot 9, Fl)R'l~^], (5) 

where Fi = F{9i) and F/ denotes the derivative to 9i. 
The flux ratios of the images follow from their magnifi- 
cations /^i, which are given by 



{PtPi 



Fa' 



Fl')R^ 



+ {a-l) a'Fl~{a^l)F.r+aF,Fl' R 



2(a-2) 



(6) 



where pi is the parity of image i. 

Since the Fourier coefficients ci and si correspond 
to a displacement of the source position (^,?7), we set 
ci = si = to remove this degeneracy. For a given 
slope < a < 2, we are left with 2n + 1 free parame- 
ters (C, ?y, C2, S2, • ■ ■ , Cn, s„) and eleven constraints (four 
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image positions and three flux ratios). For each image, 
we obtain a prediction of the source position {^i, fji) and 
magnification jli through equations ([5]) and ([S]), respec- 
tively. We find the best-fit parameters by minimizing 



where c is the speed of light and Di , Ds and Dis are the 
(angular diameter) distance to the lens galaxy, the quasar 
source and the distance from lens to source, respectively. 
The (projected) mass within the critical curve 



E 
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Mc^^cJ^ ^[a^F{9)+F" 



Rcierae, (lo) 



with A^i and Arji errors in the observed image posi- 
tions, and fj = /ii/^i+j and Afj the observed flux ratios 
and corresponding errors"'^^. During the minimization we 
guide the solution towards a smooth, realistic lens model 
by adding to x^ a constant times Em>3(cm + «„), which with Rjf" = aco/2, so that the Einstein mass becomes 



is in general close to the projected mass within the Ein- 
stein radius. The latter is the radius Re for which the 
projected mass is equal to the Einstein mass defined as 
AIe — T,cTtR^. Because Re is independent of 9, all 
higher order Fourier terms average out, and we are left 

2-a 



^m>3VHn 

suppresses (strong) deviations from an elliptic shape. 

For a range of slopes 0.5 < a < 1.5, the fits of the 
above lens model with n — A Fourier terms are shown 
in Fig. 21 Whereas the predictions of the source po- 
sition remain within the observed errors over the full 
range of slopes, the recovery of the observed radio fluxes 
places constraints on the slope. Dividing the correspond- 
ing x^ values by the ten degrees of freedom (eleven con- 
straints minus the slope parameter) provides a quality of 
fit, shown in the lower-right panel of Fig. |4l The min- 
imum value x^/10 — 1 shows that for a ~ 1.0 a good 
overall best-fit model is obtained. The corresponding 
99 %-confidence interval (dotted line) is 0.9 < a < 1.1. 

With n < 4 Fourier terms the fits are significantly 
worse (x^/10 > 2.5), and with n — 5 Fourier terms the 
fit does not improve (for n > 5 the fitting problem be- 
comes underdetermined) . If we replace the radio fluxes 
by the optical fluxes from the CASTLES website, no ac- 
ceptable lens model fit is possible (x^/10 > 30). On 
the ot her hand, the mid- infrared fluxes from lAgol et al.l 
(|2000( ) yield similar good fits as the radio fluxes. The 
best-fit slope a = 0.9 is lower because the fluxes of im- 
age A and C are respectively about 10 % and 30 % lower 
than in the radio. Image C is the faintest in both mid- 
infrared and radio, while image D is the faintest in the 
optical. Since the radio is expected to be (much) less 
affected by differential extinction and microlensing, we 
adopt the corresponding lens models. 

In the investigation below, we use the lens models with 
slopes a = {0.8,0.9,1.0,1.1,1.2}, shown in Fig. O and 
with best-fit parameters given in Table [TJ The criti- 
cal curves Rc{8) are obtained by solving the (quadratic) 
equation ^ for infinite magnification, i.e., for vanishing 
left-hand-side. The caustics then follow upon substitu- 
tion of Rc{0) in the lens equation ^. We see in Fig. [5] 
that even the two models outside the 99 %-confldence 
interval for a are relatively smooth as expected for a re- 
alistic galaxy model. 

The surface mass density of the lens models follows 
from Poisson's equation as 

S(i?, 0) = ^c\ [a^F{9) + F"ie)] R"~\ (8) 

with the critical surface mass density deflned as 

c^ D 

(9) 



47rGAA«' 

curve) equations 

at the solution i 

ward matrix inversion IIEvans fc Witt)[2003l) . 



^^ For a = 1 (flat rotation curve) equations JSj and JBJl are linear 
in the free parameters, so that the solution follows by straightfor- 



Me = EcTtI^— 1 



(11) 



Since we express all (projected) coordinates on the plane 
of the sky in arcseconds, including the radius i?, we mul- 
tiply expressions HO]) and ^^ by {Di 7r/0.648)2 to con- 
vert from arcseconds to pc for a given (angular diameter) 
distance Di to the lens galaxy in Mpc. 

From the last two columns in Table (TJ we see indeed 
that Mc « Me = 1.54 x 10^° Mq, nearly independent of 
the slope a. Taking into account an inverse scaling with 
the Hubble constant of Hq = 73kms~^ Mpc~^, our val- 
ues a re within a few per cent of previous measurements 
(e.g. iRix et al.l fl99a IWambsga nss fc Paczvnski 199j; 



Chae et al.lll998i ISchmidt et all 119981 : iTrott fc Websterl 
200i r 



4. RESULTS 
4.1. Luminous versus total mass distribution 

In the left and middle panel of Fig. [T] we show MGE 
flts (contours) to respectively the F-band and the J-band 
HST image, obtained with the software of ICappellaril 
()2002f ). while masking the quasar images. For each fit 
we assume the same position angle (PA) for the Gaus- 
sians, consistent with the adopted axisymmetry in the 
dynamical models. Although the bar and spirals cannot 
be reproduced, it provides a good description of the disk 
in the outer parts (left panel) and reproduces well the 
bulge in the inner parts (middle panel). Similarly, we de- 
rive MGE fits to the surface mass density in equation ([8]) 
of the best-fit lens models with slopes a from 0.8 to 1.2. 
In the right panel of Fig. [Tl we show the MGE fit (solid 
contours) to the surface mass density (dashed contours) 
of the overall best-fit lens model with slope a — 1.0. The 
corresponding parameters are given in Table [5] (columns 
5-7), together with the parameters of the MGE fit to 
the /-band surface brightness (columns 2-4). The latter 
parameters are "deconvolved" using a MGE model fit- 
ted to the WFPC2 PSF as determined with Tiny Tim 
(|Krist fc Hookl[200l . 

The Gaussians in the MGE fit to the /-band surface 
brightness have a PA of ~ 70°, which is similar to the 
PA of ~ 67° in the MGE fit to the surface mass density 
of lens models. B oth these va lues are consistent with the 
measurements bv lYed (|1988[ ). who found a PA of ~ 67° 



for the axis through quasar images C and D, bracketed 
by a PA of ~ 77° for the out er disk and a PA of - 39° 
for the bar (see also Fig. 1 of lTrott fc Webster! [200l ). 

As shown in Fig. [SI not only the orientation, but also 
the shape of the luminous density distribution as inferred 
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Fig. 4. — Lens model with n = 4 Fourier terms (see §[33] for details) fitted to the position of the quasar source (first two panels, in 
milliarcseconds) and to the observed radio flux ratios of the quasar images (third panel). Over the full range of slopes 0.5 < « < f .5, 
the predicted source positions from each quasar image A through D (as indicated by the different symbols) fall within the observed error 
(dotted lines). However, the recovery of the observed flux ratios (solid lines) within the observed errors (dotted lines) places constraints 
on the slope. The minimum reduced x^ value of about unity in the fourth panel shows that for a ~ 1.0 a good overall best-fit model is 
obtained, with corresponding 99 %-confidence interval (dotted line) of 0.9 ^ a < 1.1. 



1.0 

^ 0.5 

<u 
en 

D 0.0 



a = O.I 



a = 0.9 



a=1.0 



a=1.1 



= 1.2 




-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 
-Aa (arcsec) -Aa (arcsec) 



.0 -0.5 0.0 0.5 1.0 
-Aa (arcsec) 



-0.5 0.0 0.5 
-Aa (orcsec) 



-0.5 0.0 0.5 
-Aa (orcsec) 



Fig. 5. — Best-fit lens models for five different values of the slope a indicated at the top of each panel (see Fig. |4] for corresponding 
quality of the fit). The filled circles indicate the positions of the four quasar images, the cross represents the center of the lens galaxy, and 
the diamond shows the best-fit position of the quasar source. The smooth and pinched solid curves show respectively the radial critical 
curve and radial caustic. The dashed curves are the tangential counterparts when they exist (for a > 1.0). 
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Lens 


MODEL PARAMETERS 










a 


i 


V 


CO 


C2 


C3 


C4 


S2 


S'i 


S4, 


Re 


Me 


Mc 


0.8 


0.0824 


-0.0151 


2.1616 


-0.0590 


0.0001 


0.0017 


0.0593 


0.0022 


0.0030 


0.8858 


1.5414 


1.6004 


0.9 


0.0736 


-0.0143 


1.9456 


-0.0494 


0.0001 


0.0010 


0.0500 


0.0005 


0.0019 


0.8862 


1.5425 


1.5846 


1.0 


0.0651 


-0.0133 


1.7733 


-0.0415 


0.0002 


0.0006 


0.0423 


-0.0006 


0.0012 


0.8866 


1.5441 


1.5756 


1.1 


0.0570 


-0.0121 


1.6325 


-0.0349 


0.0004 


0.0003 


0.0357 


-0.0014 


0.0009 


0.8872 


1.5461 


1.5708 


1.2 


0.0492 


-0.0107 


1.5154 


-0.0293 


0.0005 


0.0002 


0.0299 


-0.0018 


0.0007 


0.8879 


1.5485 


1.5691 



Note. — Parameters of the lens models shown in Fig. [5] For each slope a, the best-fit quasar source position (5,»y) and Fourier 
coefficients Cm and Sm (ci = sq = *i = 0) are given, all in arcseconds. The last three columns provide the Einstein radius (in arcsec), the 
projected mass within this radius and within the critical curve (both in 10^" Mq). 



from the surface brightness, is very similar to that of 
the total mass density distribution as derived from the 
lens models. The left and middle panel in the top row 
show respectively the profile and the projected elliptic- 
ity e = 1 — g' of the surface mass density as function 
of the major axis on the sky-plane. The thin curves are 
for the lens models with fixed slope a (as indicated in 
the bottom-right panel). The thick solid curve is for 
the /-band surface brightness, multiplied with a constant 
mass-to-light ratio T^; = 3.4 Tq. 7. The latter is derived 
from dividing the projected mass from the lens models 
(shown in the top-right panel of Fig. [6|) within the Ein- 
stein radius, i.e., the Einstein mass Me, by the projected 



luminosity within the Einstein radius Re- 

We see that around Re (indicated by the dashed ver- 
tical line) , where the lens models are best constrained by 
the observed quasar image positions and relative flux ra- 
tios, both the slope and the ellipticity of the projected lu- 
minous and total mass density are nearly the same. The 
same holds true for the corresponding intrinsic mass den- 
sities, of which the profiles and ellipticities eintr = 1 — 9 
are shown in the left and middle panel in the bottom 
row of Fig. [6] These (analytic) deprojections (see Ap- 
pendix 1^ are under the assumption of oblate axisym- 
metry and for a given inclination which we derive from 
the flattening of the disk in the outer parts. We find 
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Fig. 6. — Radial profiles of the mass distribution of the lens galaxy in the Einstein Cross as function of radius R, in arcsec at the bottom 
and in kpc in at the top. In each panel the dashed and solid vertical lines indicate respectively the Einstein radius Re and the effective 
radius Re- The top panels shows from left to right the surface mass density profile S, the corresponding projected ellipticity e = 1 — q' , 
and the mass M enclosed within the projected radius R along the major-axis on the plane of the sky. The dotted, dash-dotted and dashed 
curves are from the best-fit lens models with (fixed) slopes a of respectively 0.9, 1.0 and 1.1 (see ii l3.3l l. The thick solid curves are from the 
observed (/-band) surface brightness, multiplied with a constant mass-to-light ratio Tg = 3.4 Tqj. The latter is obtained from dividing 
the Einstein mass (indicated by the horizontal dashed line in the top-right panel) by the (projected) luminosity within the Einstein radius. 
The bottom panels shows from left to right the mass density profile p, the corresponding intrinsic ellipticity eintr = 1 — 9i and the circular 
velocity Vdic as function of the radius R in the meridional plane. These intrinsic quantities are obtained from the MGE fits to the projected 
density, under the assumption of oblate axisymmetry and for an inclination i = 68° (see also Appendix Rt. 



TABLE 2 
Multi-Gaussian Expansion parameters 





/-band surface bri 
loglo logo-' 


jhtness 


a = 1.0 lens model 


i 


log So 


logo-' 


<?' 


1 


4.329 


-1.564 


0.700 


5.114 


-1.398 


0.645 


2 


3.935 


-0.942 


0.700 


4.556 


-1.037 


0.660 


3 


3.606 


-0.627 


0.700 


4.238 


-0.779 


0.662 


4 


3.293 


-0.239 


0.700 


3.928 


-0.544 


0.678 


5 


3.005 


-0.043 


0.700 


3.628 


-0.350 


0.673 


6 


2.845 


0.230 


0.700 


3.459 


-0.186 


0.675 


7 


2.261 


0.587 


0.700 


3.354 


-0.001 


0.675 


8 


2.160 


0.917 


0.414 


3.208 


0.220 


0.675 


9 


1.334 


1.135 


0.700 


3.069 


0.496 


0.675 


10 


- 


- 


- 


2.991 


1.000 


0.675 



Note. — The parameters of the Gaussians in the MGE fit to 
the HST/WFPC2/F814W /-band image of the surface brightness 
(columns 2-4, "deconvolved" with the WFPC2 PSF), and to the 
surface mass density of the overall best-fit lens model with slope 
a = 1.0 (columns 5-7) of the lens galaxy in the Einstein Cross (see 
also Fig. [TJ. Columns two and five give for each Gaussian com- 
ponent respectively the central surface mass density (in Mq pc~-^) 
and the central surface brightness (in Lq pc~^), columns three and 
six the dispersion (in arcsec) along the major axis, and columns 
four and seven the observed flattening. 



from the MGE fit to the F-band surface brightness (left 
panel of Fig. [T]), that Gaussian components -with a mea- 
sured flattening as small as q' ~ 0.4 are required for an 
acceptable fit. This sets a lower limit to the inclination 
of i > 66° (signi ficantly larger than i ~ 60° found by 
llrwin et al.lll989r ). Adopting a lo"wer limit for the intrin- 
sic fl attening of the disk oi q = 0.15 (e.g. iLambas et alj 
Il992f ). -we then obtain an inclination of i ~ 68°. 

Around Re the luminosity density is close to isother- 
mal (p oc R^^ intrinsic or E ex R~^ in projection) like 
the mass density of the overall best-flt lens model. To- 
wards the center the slope of the luminosity density be- 
comes shallower. The (fixed) slope of the lens model 
is not anymore (well) constrained by the quasar images 
towards the center, neither for radii a few times Re- 
At these larger radii also the luminosity density is not 
only due to the bulge, but the bar and disk start con- 
tributing, as can also be seen from the increase in the 
ellipticity. Nevertheless, within a radius R < 4", the 
intrinsic luminous and total mass distribution are very 
similar, i.e., mass follows light. This means we can in- 
fer the gravitational potential either from the observed 
surface brightness adopting a (constant) total mass-to- 
light ratio Ttot, or directly from the best-flt lens model 
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(see also Appendix El . The corresponding circular veloc- 
ity curves, defined as w^j^j. = Rd^/dR in the equatorial 
plane, are shown in the bottom-right panel of Fig. [SI us- 
ing Te = 3.4 Tq J as above. 

4.2. Axisymmetric Jeans models 

In Fig. [71 the second panel in the bottom row shows 
the square root of the combination V^ + a^ of the ob- 
served (and symmetrized) mean line-of-sight velocity V 
and dispersion a maps of the lens galaxy (Fig. [3l two 
panels on the right). The other panels in Fig. [7| show the 
line-of-sight second-order velocity moment as predicted 
by the axisymmetric Jeans models, using the solution 
given in equation (|A18p of Appendix [^ With a central 
velocity dispersio n of gp c:^ 166 km s""*^, t he Mbh ~ cr re- 
lation as given bv lTremaine et all (|2002f ) predicts a cen- 
tral black hole of mass Mbh — 6.4 x 10^ M© for the 
lens galaxy. Since the corresponding sphere-of-influence 
Mbh/c^ ~ 0.01" is not resolved by our kinematic data, 
we cannot fit for it, but do include the predicted BH 
in the Jeans models (as an additional Gaussian). Fi- 
nally, we convolve the predicted second-order velocity 
moment with the Gaussian PSF of the kinematic data 
with FWHM= 1". The (hnear) scale in kms~-^, as indi- 
cated by the color bar and limits at the top of the panel 
with the observations, is the same in all panels of Fig. [7) 

The axisymmetric Jeans models in the top row use 
the gravitational potential inferred from the /-band sur- 
face brightness, assuming that mass follows light. The 
corresponding (constant) dynamical mass-to-light ratio 
Tdyn that provides the best-fit to the observations is in- 
dicated at the top, together with the assumed value for 
the anisotropy parameter /3, increasing from left to right. 
The corresponding increase in the (constant) flattening 
of the velocity ellipsoid in the meridional plane (or de- 
crease in f^/flj), has two main effects: a less pinched 
"butterfly" shape and a weaker gradient parallel to the 
(vertically aligned) short axis. Comparing with the ob- 
servations, we see that at lower /? values the pinching is 
too strong, while at higher /3 the gradient along the short 
axis is too shallow. Therefore, it is not unexpected that 
the "best-fit" is when /3 ~ and Tdyn — 4.0 Tq ^ (in 
the /-band), but with reduced x^ — 3 significantly above 
unity. 

Part of the latter mismatch might be due to (system- 
atic) discrepancies in the data, in particular due to the 
challenging measurement of the velocity dispersion. The 
average uncertainty in the observed second-order veloc- 
ity moment is about lOkms"^, which translates into an 
error of around 12 per cent in Tdyn, or about 0.5 Tqj. 
However, the uncertainties in the observed second-order 
velocity moment vary quite a lot, from only r-^ 2kms~^ 
in the center to ~ 20kms~^ a the edge of the CMOS 
field. And although we were very careful in disentan- 
gling the contribution from the quasar and excluded the 
regions at the quasar images, the measured velocity dis- 
persion around these regions might still be systematically 
affected by (broadening due to) residual light from the 
quasar source. This possibly explains the observed excess 
in the observed second-order velocity moment around the 
quasar images with respect to the accurately measured 
central value. At the same time, such an excess might 
mimic a steeper gradient in the observed second-order 



velocity moment. Indeed, in all Jeans models in the top 
row of Fig. [7| the prediction at the center is significantly 
higher than the observed value. As a result, we expect 
that the dynamical mass-to-light ratio is overestimated. 

If, alternatively, we infer the gravitational potential 
directly, i.e., without the need for a mass-to-light ratio, 
from the (deprojected) surface mass density of the best- 
fit lens models, the axisymmetric Jeans models predict 
maps of the second-order velocity moments that are sim- 
ilar in shape but with lower values near the center. This 
is shown in the middle row of Fig. [7] for the overall best- 
fit lens model with slope a = 1, and for the same range 
in anisotropy parameters /3. The formal best-fit is ob- 
tained for /3 ~ 0, but the observed central value is more 
closely matched in case of mild velocity anisotropy in the 
meridional plane with /3 > 0.1. For a = 0.9. i.e., den- 
sity slopes steeper then isothermal, the observed central 
value is best matched when /? ~ 0, with at the same 
time a steeper gradient along the short-axis, as shown 
in the first panel in the bottom row of Fig. [71 The last 
two panels show that for a > 1, velocity isotropy in the 
meridional plane (fi — 0) results in a central value much 
lower than the observed second-order velocity moment. 

We conclude that for the lens models with slopes 
a = 1.0 ± 0.1, as constrained by the lensing geometry, 
and for a mild anisotropic velocity distribution in the 
meridional plane with /? = 0.1 ± 0.1, the predictions of 
the second-order velocity moments match the accurately 
measured central value, and follow the gradient within 
the observational uncertainties outside the region that 
might be affected by (residual) quasar light. The lat- 
ter results in broadening of the velocity dispersion, so 
that the best-fit dynamical mass-to-light ratios we derive 
when we infer the gravitational potential from the surface 
brightness, are likely overestimated. Fitting instead only 
the unaffected and accurately measured second-order ve- 
locity moment in the center, we find for /3 = 0.1 a best- fit 
value of Tdyn — 3.7 Tq/, with an error around 0.5 Tqj. 
On the other hand, both the Einstein mass and luminos- 
ity within the Einstein radius are well constrained, re- 
sulting in at most a few per cent error in T^ — 3.4 Tq^. 
Therefore, in what follows, we adopt the latter as the 
total mass-to-light ratio in the inner bulge-dominated re- 
gion of the lens galaxy. 

In the above axisymmetric Jeans models we fix the 
inclination to z = 68° determined from the observed el- 
lipticity of the outer disk (see § 14. 1|) . When we vary the 
inclination over a reasonable range above the lower limit 
i > 66°, there is no significant change in the above re- 
sults, which is expected because the inclination only has 
a weak effect on the ma. ss-to-light ratio (see also Fig. 4 
of lCappellari et al.l[200^ . 

4.3. Stellar mass-to-light ratio 

We have shown that in the bulge-dominated inner re- 
gion of the lens galaxy, the total mass distribution closely 
follows the light distribution. This does not mean there 
is no dark matter present, since if dark matter follows 
(nearly) the same distribution it might contribute to part 
of the above estimated total mass-to-light ratio Ttot — 
3.4 Tq J. To determine this possible constant dark mat- 
ter fraction, we need to measure the stellar mass-to-light 
ratio T*. To this end we use the single ste llar popu- 
lation (SSP) models of IVazdekis et al.l ()1996D . adopting 
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Fig. 7. — Axisymmetric Jeans model of the lens galaxy in the Einstein Cross. The second panel in the bottom row shows the square 
root of the combination V^ +a^ of the observed mean line-of-sight velocity V and dispersion a maps (Fig.|3]l. The (linear) scale in kms~^ 
is indicated by the color bar and limits at the top, and is t he sa me for all oth er pa nels, which show the line-of-sight second-order velocity 
moment as predicted by the axisymmetric Jeans models (§[32] and Appendix IA.3| l. In the top row, the gravitational potential is inferred 
from the /-band surface brightness, assuming that mass follows light. The corresponding (constant) dynamical mass-to-light ratio T^yn 
that provides the best-fit to the observations is indicated at the top, together with the assumed value for the anisotropy parameter /3, 
increasing from left to right. The latter parameter represents the (constant) flattening of the velocity ellipsoid in the meridional plane. 
The middle row follows the same sequence of increasing /3 values, but the gravitational potential is now computed directly, i.e., without 
mass-to-light ratio conversion, from the overall best-fit lens model with slope a = 1.0. In the bottom row, besides the observations in 
the second panel, the predictions for different slopes a are shown, with /3 = 0, corresponding to an isotropic velocity distribution in the 
meridional plane. 

Next, we convert the resulting magnitudes in the F555W, 
F675W, F814W (WFPC2), F160W and F205W (NIC- 
M0S2) filters to (rest frame) magnitudes in respectively 
the Johnson V , J, R . H and K filters. \ ye use the method 
described in § 3 of Ivan de Ven et all ()2003f ) to take si- 
multaneously into account the filter responses and the 
redshift of the lens galaxy. 

The resulting colors are best fitted by a SSP model 
with age t ~ 8Gyr and metaUicity [Fe/H] ~ 0.2, corre- 
sponding to an /-band stellar mass-to-light ratio T^, ~ 
3.3 Tq J. The lu confidence limits yield a range in t from 
7 to 14Gyr, and in [Fe/H] from 0.0 to 0.3, corresponding 
to a range in T^, from 2.8 to 4.1 Tq_/. The lower limit im- 
plies at most ~ 20 per cent of dark matter within i? < 4", 



as initial reference the 'Kro upal ()2001[ ) initial mass func- 
tion (IMF) with lower mass cut-off 0.01 Mq (faintest star 
0.09 Mq) and upper mass cut-off of 120 Mq. 

We derive a high-S/N spectrum from our data cube by 
collapsing the central spectra that are unaffected by the 
quasar images. Even so, the Ca II triplet is equally well 
fitted by SSP models of nearly all ages and metallicities 
of about [Fe/H] = ±0.2 around solar, which leaves the 
stellar mass-to-light ratio nearly unconstrained. There- 
fore, we also derived colors from archive HST images 
(PI: Kochanek). After matching all images, using both 
the quasar positions and isophotes of the lens galaxy, we 
compute the magnitudes in a circular aperture of radius 
R = 0.2" (well within the quasar images at i? ~ 0.9"). 



12 



van de Ven et al. 



but the results are fully c onsistent with no dark matter 
at all. The SSP models of lMarastonI (|1998l I2005D with a 
Kroupa IMF yield similar results. Adopting a iSalpeten 
([1955) IMF with the same lower and upper mass cut-offs, 
increases T^, by about 30 per cent. This not only implies 
no dark matter, but even for the lower limit T^, > Ttot, 
which is unphysical. Th is is consistent with ev idences for 
both late-type galaxies ('Bell fc de J onEJr2001[ ) and early- 
type galaxies (Cappellari et al. 2006) that, if the IMF is 
universal, a Salpeter IMF is exclud ed, whereas a Kroupa 
IMF matches the observations (see lde Jong fc Belll2007l 
for a review). 

5. DISCUSSION AND CONCLUSIONS 

We used the GMDS-North integral-field spectrograph 
to obtain two-dimensional stellar kinematics of the lens 
galaxy in the Einstein Cross. In addition to the four 
bright quasar images and the distance of the lens galaxy 
{Di = 155 Mpc), in particular the presence of sky lines 
in the observed Ca II triplet region made the extraction 
of the absorption line kinematics challenging. Even so, 
we were able to derive high-quality line-of-sight velocity 
V and dispersion a maps of the bulge-dominated inner 
region i? < 4", reaching about two-thirds of the effective 
radius Re — 6" of this early-type spiral galaxy. The V 
map shows regular rotation up to ~ 100 km s^^ around 
the minor axis of the bulge, consistent with axisymmetry. 
The (T map shows a weak gradient increasing towards a 
central value of 166 ± 2kms^^. 

The only other direct (single) measure ment of the 
veloci ty dispersion is 215 ± SOkms"'^ by iFoltz et aTl 
(|1992t ). For a singular isothermal sphere lens model, 
we can use the relati on A9 = 8n{asis/cy Dis / Dg (e.g. 
iKochanek et al.l l2000f ) with a separation A6 = 1.8" of 
the four quasar images, to obtain a simple estimate for 
the dispersion of ctsis — 180kms~^. Taking into account 
aperture correction and a range in velocity anisotropy, 
Ivan de Ven et al.l ()2003f ) convert this to a central veloc- 
ity dispersion of 1 68 ib 17kms~^. The K ing and de Vau- 
couleurs models of lKent fc F alco' (1988") predict a similar 
value of ~ 166kms~^, and also Barnes et al. (1999) find 
a value of 165±23kms~^ based on their two H I rotation 
curve measu rements. It i s like ly that the long-slit mea- 
surement by iFoltz et all (|1992t l is affected by the bright 
quasar images, whereas our spatially resolved measure- 
ments allow for a clean separation of the quasar contri- 
bution. 

A large variety of different lens models have been con- 
structed for the Einstein Cross, most of which fit the po- 
sitions of the quasar images but not their relative flux ra- 
tios. Adopting the scale-free lens model of Evans fc WittI 
(J2003i) . we found that fitting at the same time also the 
(radio) flux ratios constrained the slope of the total mass 
surface density E oc R"~^ to be a = 1.0 ± 0.1. The to- 
tal mass within the Einstein radius Re = 0.89", i.e., the 
Einstein mass, is Me = 1-54 x 10^" M©, nearly indepen- 
dent of the slope a, and consistent with previous mea- 
surements in the literature. Dividing by the projected 
luminosity within Re, as measured from the observed 
/-band surface brightness, we obtained a mass-to-light 
ratio of T^; = 3.4 Tqj, with an error of at most a few 
per cent. 

We determined the mass-to-light ratio in an addi- 
tional, independent way by fitting dynamical models to 



the (combined) observed V and a maps. We used the 
solution of the axisymmetric Jeans equations to pre- 
dict this second-order velocity moment, with the grav- 
itational potential inferred from the deprojected /-band 
surface brightness, assuming that mass follows light. We 
expect the resulting best-fit constant mass-to-light ratio 
to be an overestimation due to possible residual contribu- 
tion from the quasar light. Even so, restricting the fit to 
the unaffected and accurately measured central second- 
order velocity moment, we found Tdyn = 3.7 ± 0.5 Tqj. 
When we used instead the gravitational potential that 
follows directly, i.e., without mass- to- light ratio conver- 
sion, from the surface mass density of the best-fit lens 
models, we arrived at a similar prediction of the second- 
order velocity moment. We showed that the reason is 
that the luminous and total mass distribution, as inferred 
from respectively the surface brightness and lens model, 
are very similar. This implies that the inner region of 
the lens galaxy mass (closely) follows light, with a total 
mass-to-light ratio Ttot — 3.4 Tqj. 

By fitting single stellar population models to measured 
colors of the center of the lens galaxy, we estimated an /- 
band stellar mass-to-light ratio T^, from 2.8 to A.ITqj. 
Although a constant dark matter fraction of 20 per cent 
is thus not excluded, it is likely that dark matter does 
not play a significant role in the inner region of this 
early-type spiral galaxy of luminosity ^ L-k. This is 
consistent with indications that loss-lumi nous early-type 
galaxies (e.g. [Ge rhard et al. 2001; Cappel lari et al.ll2006l : 
iKoopmans et al. 2006; Thomas ct al. 2007) and lat e-type 
galax i es (e.g. Paluiias fc Williams 2000; Bell fc d e Jond 
120011: iKassin et all I2006D lack the dark matter that is 
ubiquitous in dwarf galaxies as well as in (at least the 
outer parts of) giant elliptical and spiral galaxies. 

The constraint a = 1.0 ± 0.1 on the slope of the lens 
model implies that the intrinsic total mass density is close 
to isotherma l, consistent with prev ious studies of lens 
galaxies fe.g. iKoopmans et al.ll20(56r ). However, one has 
to be very careful not to over- interpret this result, not 
only due to assumptions in the modeling, but most of 
all because the slope is only (well) constrained around 
the Einstein radius Re- If indeed mass follows light in 
the inner region of this lens galaxy, the (deprojected) 
surface brightness indicates deviations from isothermal 
towards the center where the slope becomes shallower as 
well as a possible steepening at larger radii. Interest- 
ingly, the construction of realistic galaxy density profiles 
with a stellar and dark matter component that are both 
non-isothermal, shows that the combined slope and cor- 
responding lensing properties a re nevertheless consistent 
with isothermal around Re ()van de Ven et al.l l2008bl : 
iMandelbaum et al]|2008D . With more extended images, 
e.g. in the case of galaxy-galaxy lensing, one can place a 
stronger constr aint on the total mass distribution of the 
lens g alaxy (e.g. lBarnabe fc Koopman j|2007t Isiivu et al.l 
|2008[) . Since the extension in these cases is still mostly 
tangential, constraining a large radial range is in partic- 
ular possible when lensing occurs in multiple (redshift) 
planes, resulting in images at different Einstein radii (e.g. 
iGavazzi et al]|2008[ ). 

Nevertheless, already the quasar images provide an ac- 
curate constraint on the total mass within Re, nearly in- 
dependent of the details of the lens model (e.g. KochaneM 
Il991t lEvans fc WittI [200l . Therefore, at least around 



Stellar dynamics and gravitational lensing 



13 



this radius no higher-order velocity moments are needed 
to break the mass-anisotropy degeneracy. Towards the 
center the surface brightness increases steeply, so that 
high enough S/N to measure velocity moments beyond 
V and a might be achievable even at higher redshift. In 
the outer parts, the degeneracy might be (partially) bro- 
ken by combining the kinematics of stars and/or discrete 
tracers, such as globular clusters and planetary nebulae, 
with total mass estimates fror n. hot X-ray gas an d/or 
weak lensing measurement (e.g. iGavazzi et al]|2007l ). In 
the current and even more in the upcoming extensive 
and deep photometric surveys, numerous (strong) gravi- 
tational lensing systems will be discovered. They provide 
important and independent constraints on the total mass 
distribution in galaxies, especially in combination with 
(resolved) st ellar kinematics, as we show e d in this pa- 
per ( see also iBarnabe fc KoopmansI 120071 : ICzoske et al.l 
[200l . 
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APPENDIX 
MULTI-GAUSSIAN EXPANSION 



We summarize the (numerically) convenient expressions of the axisymmetric intrinsic density and gravit ational 
poten tial in the case of a Multi-Gaussian Expansion (MCE) of the surface density (Monnct et al. 1992; EmscUcm et al.l 
1 19941) . Next, we show that also the gravitational lensing properties can be readily computed, so that the MGE method 
can be used to efficiently construct general lens models. Finally, we give the line- of-sight second-ord er velocity moment 
as a solution of the axisymmetric Jeans equations which is discussed in detail in iCappellaril ()2008[ ) . 



Axisymmetric density and potential 
We parameterize the surface brightness I{x',y') by a sum of Gaussian components 
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each with three parameters: the central surface brightness /qj, the dispersion cr' along the major z'-axis and the 
flattening g'. In case of an (oblate) axisymmetric system viewed at an inclination i > 0, the corresponding intrinsic 
luminosity density is 
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with intrinsic dispersion aj = tr' and intrinsic flattening qj given by q'j sin^ i = q 

The intrinsic mass density follows as pj = Tj Vj , where the mass-to-light ratio Tj per Gaussian component is a 
free parameter. Alternatively, if the surface mass density Yi{x' ,y') is available (e.g. from gravitational lensing), an 
MGE-parametrization as in (JA1|) directly yields pj as in (|A2|1 . but with Iqj in both expressions replaced by the central 
surface mass density Eo J . 

The corresponding gravitational potential follows upon (numerical) evaluation of (|Emsellem et al.lll994 eq. 39) 
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where we have introduced 



Tj{u) = exp<^ -— 2 
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and Q^ (u) = 1 — (1 — g^)^^. The total mass per Gaussian component is given by Mj 



27ro-' g'Soj, with So,, = T,- /q, 



when the mass density is inferred from the surface brightness. 

In the latter case, we have implicitly assumed a constant mass-to-light ratio T j per Gaussian component by taking 
it outside the integral. Nevertheless, we can still mimic a (radially) vary ing mass-to-light ratio by considering the Tj 
of the Gaussian components as free parameters (see e.g. Ivan de Ven etal . 2006; van den Bosch et al. 2006). However, 
it is common in dynamical studies of the inner parts of galaxies (e.g. ICappellari et al.ll2006) to assume the total mass- 
to-light ratio to be constant, i.e., T j — Ttot for each Gaussian component j. Since Ttot may be larger than the stellar 
mass-to-light ratio T^,, this still allows for possible dark matter contribution, but with a constant fraction. 
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Axisymmetric lens model 

Under the thin-lens approximation, the gravitational lensing properties of a galaxy are characterized by the deflection 
potential 4'{x' , y') and its (partial) derivatives. The deflection potential follows from projecting the potential along the 
line-of-sight or by solving the two-dimensional Poisson equation V^0 = 2k. Here, k = S/Sc is the normalized surface 
mass density with the critical lensing value Sc defined in equation ©. 

Similar to the surface brightness in equation (jAip . we parameterize this so-called "convergence" k{x' ,y') by a sum 
of Gaussian components 
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The corresponding deflection potential is then 
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where Mk = Sttct^ q^So,fc is the total mass per Gaussian component, and we have introduced 
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andQ',2(w) = l-(l-gf)«2^ 
The first-order partial derivatives of 0fc(a;', y') follow as 
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whereas the second-order partial derivatives are given by 
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The above single integrals can be readily evaluated numerically, so that the lensing properties follow in a straightforward 
way. The image positions {x',y') are related to the source position (^,77) by the lens equation 
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Given the parity p of the image, the magnification /i is given by 
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The Einstein radius Re follows from solving k{Re) = 1, where k{R') is the average convergence within the projected 
radius R'. After substituting x' — R' cos 6' and y' — R' snvB' in equation (jASp and performing the integral over i?', we 
are left with 
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where V'iO') = cos 6' + sin 9' /q'l^ . The mass within the Einstein radius then follows as AIe = Sc7ri?|j. 



Axisymmetric Jeans equations 

We consider a luminous component with Gaussian intrinsic luminosity density ja,- given by (IA2p . of which the 
observed kinematics trace the unde rlying gravitational potential $ — J2k^k, with the sum over all (luminous and 
dark) components ^k given by (jA3p . 
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Assuming that the velocity ellipsoid is aligned with the cylin drical {R, (f>, z) coordinate system, so that vrVz 
can readily solve Jeans equation ^ as (jEmsellem et al.lfi994 eq. 42) 
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with Tk{u) defined in (jA4[) and 
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= 0, we 
(A13) 
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If we next also assume a constant flattening of the velocity ellipsoid in the meridional plane, we can write [f^]j = &j [v1]j, 
and solve Jeans equation ([!]) as 
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where we have introduced 



GM^^) = 1 - b,qj - Ml - qi) + (1 - b^j)F,ky. (A16) 

Since the constant bj may be different for each luminous component j, the total anisotropy in the meridional plane 
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is allowed to vary throughout the system. I n case of isotropy in the meridional plane all bj = 1 and thus /? = 0, 
corresponding to a two-integral DF f{E, Lz) (jEmsellem et al.lll994[ ). 

Aft er substitution of these expressions in equation ([3]), the line-of-sight integral can be solved, resulting in (|Cappellaril 
[20081 eq. 27) 
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where A and B are functions of u defined as 
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The remaining single integral can be readily evaluated numerically. 

This provides a direct prediction of the combination V^ + a^ of the observed mean line-of-sight velocity V and 
dispersion cr at a given position {x',y') on the sky-plane, through the (luminosity weighted) sum v^^^{x',y') = 
Sj -(jhoslj/ Si ^j over the corresponding luminous components. 
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